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SUMMARY  AND  CONCLUSIONS 


An  accurate  full-wave  method  for  determining  the  effect 
on  ELF  radio  propagation  of  a localized,  cylindrically- 
symmetric  disturbance  has  been  developed.  Application  of 
this  method  to  the  disturbance  produced  by  a high-altitude 
nuclear  explosion  gives,  when  the  transmitter  is  outside  of 
the  disturbed  region,  greater  attenuation  in  most  of  this 
region  than  does  the  two-dimensional  WKB  approximation  used 
by  previous  workers,  but  less  attenuation  at  the  far  edge  of 
the  region.  For  example,  10  min  after  a typical  megaton-range, 
high-altitude  nuclear  explosion,  the  signal  level  at  the  far 
edge  is  calculated  to  be  3 dB  below  that  before  the  bui  t, 
while  the  two-dimensional  approximation  gives  a 7-dB  decrease. 
Such  differences  are  very  important  for  ELF  systems  whose 
enormous  power  requirements  may  necessitate  operation  at 
minimal  signal-to-noise  ratios. 

The  analysis  in  this  report  also  shows  that  ELF  fields 
at  the  ground  are  sensitive  to  the  details  of  ionospheric 
conductivity  profiles  only  in  two  rather  limited  altitude 
ranges.  The  lower  of  these  is  the  altitude  range  where  the 
conduction  and  displacement  currents  are  comparable.  The 
upper  region  is  where  the  conductivity  scale  heioht  is 
comparable  with  the  skin  depth.  Future  effort  in  obtaining 
improved  conductivity  profiles  for  ELF  propagation  calcula- 
tions should  therefore  focus  on  these  two  altitude  ranges. 
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SECTION  I 


INTRODUCTION 


The  purpose  of  this  report  is  to  present  a method  for 
determining  the  effect  on  ELF  (extremely  low  frequency)  radio 
communication  of  a localized,  cylindrical ly-symmetric  iono- 
spheric disturbance,  such  as  would  be  produced  by  a solar 
proton  event  or  a high-altitude  nuclear  burst.  At  ELF, 
wavelengths  are  frequently  comparable  with  the  dimensions  of 
ionospheric  disturbances,  either  natural  or  artificial. 
Consequently,  approximate  methods  for  handling  such  problems 
at  VLF , such  as  the  WKB  approximation  or  the  Fresnel  diffrac- 
tion model  of  Crombie  [1],  are  not  applicable  at  ELF.  At 
these  lower  frequencies,  it  becomes  necessary  to  work  with 
full-wave  solutions  of  Maxwell's  equations  in  the  earth- 
ionosphere  waveguide. 

There  are  well-known  techniques  for  obtaining  full-wave 
solutions  for  an  earth-ionosphere  waveguide  in  which  there  is 
inhomogeneity  in  only  the  vertical  direction  [2-4].  For 
certain  idealized  conductivity  profiles,  analytic  solutions 
are  possible.  However,  for  more  realistic  profiles,  solutions 
can  be  obtained  only  by  numerical  methods.  Either  procedure 
involves  the  solution  of  an  eigenvalue  equation,  each  eigen- 
value and  corresponding  eigenfunction  corresponding  to  a 
particular  waveguide  mode.  The  electromagnetic  field  radiated 
by  any  source  in  the  waveguide  can  be  expressed  as  a super- 
position of  waveguide  modes,  with  relative  amplitudes  depending 
on  the  nature  of  the  source.  At  ELF,  this  representation  of 
the  electromagnetic  field  is  especially  convenient  because  only 
the  lowest,  or  TEM,  mode  is  nonevanescent  in  the  earth-ionosphere 
waveguide.  At  distances  from  a source  exceeding  2 or  3 times 
the  height  of  the  waveguide,  only  its  TEM  component  remains. 

The  source  of  particular  interest  here  is  a horizontal 
electric  dipole  antenna  located  at  or  near  the  ground.  For  a 
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horizontal  dipole  with  current  moment  Id£,  oriented  along  the 
direction  a = 0,  and  varying  with  time  as  e la)t,  the  vector 
potential  for  the  TEM  mode  is 

U IdJi  . 

Az  = -fjj—  Ac  cosa  f0(z)H}i)  (kSQr)  (1-1) 


where  k = io/c  is  the  wave  number,  n is  the  index  of  refrac- 

g 

tion  of  the  ground,  z is  the  height  above  the  ground,  and  r 
is  the  horizontal  distance  from  the  source.  The  quantity  SQ 
is  an  eigenvalue  characterizing  the  TEM  mode  and  fQ(z)  is  the 
corresponding  eigenfunction  normalized  to  fQ  = 1 at  the  ground 
("height-gain  function").  The  quantity  A , which  has  the 
dimensions  of  a reciprocal  length,  is  known  as  the  "excitation 
factor. " 

Sufficiently  far  from  the  source,  the  Hankel  function  in 
Equation  (1-1)  may  be  approximated  by  its  large  argument 
expansion,  which  has  the  form 


H 


(1) 

1 


(kSor) 


1 

r1/2 


ikS  r 
o 
e 


(1-2) 


From  this  it  can  be  seen  that  the  real  part  of  the  (complex) 
eigenvalue  SQ  is  the  reciprocal  of  the  horizontal  phase  velocity 
in  units  of  c.  The  imaginary  part  of  SQ  is  related  to  the 
horizontal  attenuation  rate  by 

a = 8.7  x 103  k!m(S  ) dB/(1000  km)  (1-3) 


where  k is  the  wave  number  in  km 

When  the  properties  of  the  waveguide  are  laterally,  as 
well  as  vertically,  inhomogeneous,  such  as  in  the  presence 
of  a localized  disturbance,  the  vector  potential  no  longer  has 
the  simple  form  of  Equation  (1-1) . However,  if  the  scale 
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lengths  for  the  lateral  gradients  are  much  larger  than  those 
for  the  vertical  gradients,  which  is  almost  always  the  case, 
it  can  be  expected  that  the  vertical  variation  of  the  vector 
potential  will  be  determined  primarily  by  the  local  ionosphere. 
The  lateral  variation  of  the  vector  potential  at  the  ground  is 
then  governed  by  a two-dimensional  Schroedinger  wave  equation 
in  which  the  effects  of  the  localized  disturbance  appear  in 
the  form  of  a complex  localized  potential.  The  problem  then 
reduces  to  the  determination  of  the  scattering  from  such  a 
potential.  These  concepts  form  the  basis  of  the  theory  which 
will  be  developed  in  the  following  sections. 

A related  approach  to  the  problem  has  been  developed  by 
lcox  [5].  He  presents  a descriptive  semiquantitative  method 
partial  wave  scattering  from  a cylindrical ly-symmetric 
potential  from  which  upper  bounds  are  obtainable  on  the 
absolute  magnitude  of  scattering  and  absorption  of  ELF  waves 
in  a horizontally  perturbed  earth-ionosphere  waveguide.  The 
method  presented  here  allows  a calculation  of  the  electro- 
magnetic field  everywhere  in  such  a waveguide. 
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SECTION  II 


r 


EFFECT  OF  LATERAL  GRADIENTS  ON  TEM  MODE 

The  equations  which  govern  electromagnetic  propagation 
in  the  earth-ionosphere  waveguide  are  Maxwell's  equations  for 
a conducting  medium  together  with  a generalized  Ohm's  law 

1 = OE  + js  (2-1) 

where  ;js  is  the  source  current  density.  Due  to  the  presence 
of  the  earth's  magnetic  field,  the  medium  is  anisotropic  and 
the  conductivity  o is  in  general  a tensor.  However,  below 
about  70  km  magnetic  effects  become  unimportant  and  the  medium 
is  essentially  isotropic.  The  conductivity  below  this 
altitude  can  therefore  be  considered  a scalar.  It  will  be 
shown  later  that  the  TEM  field  radiated  by  a source  in  the 
waveguide  falls  off  very  rapidly  above  the  altitude  at  which 
the  skin  depth  becomes  comparable  with  the  local  conductivity 
scale  height.  In  a nuclear  environment,  this  altitude  is 
generally  below  70  km.  For  ambient  ionospheric  conditions, 
this  altitude  is  near  70  km  in  the  daytime  but  above  70  km 
at  night.  In  the  following,  it  will  be  assumed  that  the 
ambient  ionosphere  is  a daytime  one,  so  that  magnetic  effects 
may  be  neglected  in  the  undisturbed  as  well  as  the  disturbed 
region  of  the  ionosphere. 

The  case  to  be  considered  is  that  of  a localized  cylin- 
drically  symmetric  disturbance  of  an  otherwise  laterally 
homogeneous  daytime  ionosphere.  It  is  convenient  to  use  a 
cylindrical  coordinate  system  (p,4>,z)  about  the  axis  of  the 
disturbance,  with  the  origin  at  the  ground.  With  the  above 
assumptions,  and  an  assumed  time  dependence  e lu)t,  Maxwell's 
equations  become 


6 


x E = iu)B 


->■  -> 

V x B 


yQ [a ( p , z) 


i^lE  + P0DS 


(2-2) 


(2-3) 


The  second  of  these  may  be  written 


V x B = - An2,(.P,zIk!  g + y - 
U)  oJs 


where  k = u)/c  is  the  free-space  wave  number  and 


(2-4) 


2 . . i . io ( p , z) 

n ( p , z)  = 1 + - K ' — 


(2-5) 


is  the  square  of  the  local  index  of  refraction. 

It  is  convenient  to  express  the  fields  in  terms  of  a 

— ► 

vector  potential  A and  a scalar  potential  ip . Thus  one  writes, 
in  the  usual  manner, 


B = $ x A 


E = + iu>A 


(2-6) 


(2-7) 


Substitution  of  Equation  (2-6)  and  (2-7)  into  Equation  (2-4) 
leads  to 


V x 7 x A = (V(V’A)  - V2  A]  = 2—  k2($ip  + icoA)  + u 3 


If  the  gauge  of  A is  chosen  to  satisfy  the  condition 


lui  3 + . 2 . 

— ^ V • A = k ip 

n 


(2-8) 


(2-9) 


the  preceding  equation  becomes 


2 -*  2 1 2 2* 

V A + n i v • A)  V ^2  + n * a = 

n 


" Mo3s 


(2-10) 


For  a horizontal  dipole  source,  A has  components  in  the  direc- 
tion of  j and  in  the  z-direction.  However,  if  lateral  gradients 
2 s 

of  n are  neglected  compared  with  vertical  gradients,  only  the 
■> 

z-component  of  A remains  finite  at  distances  from  the  source 
exceeding  the  thickness  of  the  waveguide.  Thus,  for  such 
distances,  Equation  (2-10)  with  appropriate  boundary  conditions 
can  be  satisfied  by  a vector  potential  of  the  form 


A = (0,0, Az) 


(2-11) 


The  wave  equation  then  reduces  to 


a 1 z\  1 2 A . .2. 

t: — — or  t: + — v . A 4-  k A 

d z L2  3z  2 i z z 

\ n / n 


(2-12) 


where 


n2  _ 1 S / 3 \ 1 3 

V1  D 3p  P3p)  2 2 

P 9$ 


(2-13) 


is  the  transverse  Laplacian. 

2 

In  the  case  of  a laterally  homogeneous  medium,  where  n is 
a function  only  of  z,  the  vertical  and  transverse  parts  of 
Equation  (2-12)  can  be  separated.  Thus,  one  writes 


A ( p , i}> , z ) = g(p,<J')f(z] 


(2-14) 


which  leads  to  the  two  eigenvalue  equations 
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0 


(Bessel's  equation) 


(2-15) 


2 2 2 
Vfg+  k^g  = 


d_  !l_  df 
dz  I 2 dz 
\n 


+ k‘ 


(2-16) 


where  S is  a constant.  The  various  TM  modes,  of  which  the  TEM 
is  the  lowest,  are  solutions  of  this  form  satisfying  appropriate 
boundary  conditions.  The  only  constraint  on  the  solutions  of 
Equation  (2-15)  is  the  radiation  condition,  which  requires 
outgoing  waves  at  large  lateral  distances.  Such  solutions 
exist  for  any  value  of  the  parameter  S.  However,  the  function 
f(z)  must  satisfy  two  boundary  conditions,  one  at  the  ground 
and  the  other  at  large  altitude.  These  can  be  satisfied 
simultaneously  only  for  certain,  in  general  complex,  values 
of  S.  The  eigenfunctions  corresponding  to  these  eigenvalues 
are  the  TM  modes  of  the  waveguide. 

In  the  presence  of  lateral  gradients,  a separation  of 
Equation  (2-12)  in  the  form  of  Equation  (2-14)  is  not  possible. 
However,  since  the  scale  lengths  for  the  horizontal  variations 
of  the  ionospheric  index  of  refraction  are  a few  orders  of 
magnitude  larger  than  those  for  vertical  variations,  it  seems 
reasonable  to  assume  that  the  vertical  dependence  of  the  vector 
potential  is  determined  primarily  by  the  local  ionosphere. 

Thus,  it  is  assumed  by  analogy  with  Equation  (2-16)  that  the 
vertical  variation  of  Az  is  governed  by  the  equation 


a 

3z 


9A 


2 , . az 

i n (p,z) 


+ k‘ 


1 - 


S2(p) 
n2 ( p , z) 


A = 0 
z 


(2-17) 


2 

where  S (p)  is  the  eigenvalue  corresponding  to  a laterally 
homogeneous  waveguide  with  the  vertical  conductivity  profile 
that  exists  at  p. 
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Near  the  perfectly  conducting  ground,  Az  has  the  form 

2 

Az  = g(p»<t>)  + G < p ,4> ) + ...  (2-18) 

2 

Since  n = 1 near  the  ground,  Eqs.  (2-12),  (2-17)  and  (2-18) 

give 

G ( p , <j> ) = - k2  (l-S2  ( p))  g ( p , 0 ) (2-19) 

V2  g ( p , <t> ) + k2S2  ( p)  g(p,<j>)  = 0 (2-20) 

The  solution  of  the  problem  requires  first  a determination 

2 

of  S (p) , by  a solution  of  Equation  (2-17)  at  as  many  lateral 

distances  as  necessary,  followed  by  an  integration  of 

Equation  (2-20) . The  first  step  can  of  course  by  carried  out 

by  one  of  the  techniques  mentioned  in  the  introduction  for 

the  numerical  solution  of  an  eigenvalue  equation.  However, 

this  could  involve  a large  amount  of  computation,  which  is 

probably  unwarranted  in  view  of  the  considerable  uncertainties 

in  ionospheric  conductivity  profiles.  As  an  alternative,  an 

2 

approximate  method  for  determining  S (p),  involving  very 
little  calculation,  will  be  presented  in  the  next  section. 

It  will  be  shown  that  approximate  eigenvalues  obtained  in 
this  manner  agree  very  well  with  those  calculated  by  full  wave 
solutions,  and  are  therefore  quite  adequate  for  the  present 
purposes . 
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SECTION  III 

APPROXIMATE  METHOD  FOR  EIGENVALUE  DETERMINATION 

The  method  for  determining  approximate  eigenvalues  depends 
on  the  fact,  to  be  demonstrated  presently,  that  the  eigenvalue 
depends  on  the  details  of  the  conductivity  profile  only  in 
two  limited  altitude  ranges.  The  lower  of  these  is  the  neigh- 
borhood of  the  altitude  h at  which  a = e w.  The  upper  reqion 

° ° 2-1 
is  the  neighborhood  of  the  altitude  h^  at  which  a = (4you)C^)  > 

where  is  the  conductivity  scale  height  at  the  altitude  h^. 

The  method  consists  first  of  approximating  the  conductivity 

in  each  region  by  an  exponential  with  the  local  conductivity 

scale  height.  Two  analytic  solutions  of  Maxwell's  equations 

are  then  obtained.  One  obeys  the  proper  boundary  condition 

as  z -*■  00  and  is  valid  in  the  altitude  range  hQ  <<  z < °°.  The 

other  obeys  the  appropriate  boundary  condition  at  the  ground 

and  is  valid  in  the  altitude  range  0 < z <<  h-^.  There  is  a 

region  hQ  <<  z <<  h^  where  both  solutions  are  valid.  The 

eigenvalue  S(p),  which  appears  as  a parameter  in  these 

solutions,  is  determined  by  requiring  that  the  solutions 

agree  in  the  overlap  region. 

The  basic  equations  are  Equation  (2-9),  which  becomes 


iu) 

2 

n 


k2^ 


(3-1) 


and  Equation  (2-17) . These  equations  must  be  solved  subject 
to  appropriate  boundary  conditions  at  the  ground  and  at  large 
altitude.  At  ELF,  the  ground  can  be  considered  perfectly 
conducting,  and  the  horizontal  electric  field  at  the  ground 
must  therefore  vanish.  At  very  large  altitudes,  the  solution 
is  required  by  the  radiation  condition  to  contain  only  upgoing 
waves . 
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The  conductivity  a is  an  increasing  function  of  altitude. 
The  altitude  hQ  at  which  the  conduction  current  becomes  equal 
to  the  displacement  current,  i.e.,  at  which  a(hQ)  = eQu), 
roughly  marks  the  lower  boundary  of  the  ionosphere.  For 
altitudes  z >>  hQ, 

n2  = — >>  1 (3-2) 

eo“ 


and  Equation  (2-17)  becomes  approximately 


„ e a)  9A  \ 0 

f-  + = 0 

3z  \io  3z  / z 


(3-3) 


Combined  with  Equation  (3-1),  this  leads  to  the  wave  equation 


— | + ip  au)ip  = 0 
dz^  ° 


(3-4) 


In  any  limited  altitude  range,  the  conductivity  may  be 
approximated  by  an  exponential,  with  a scale  height  appro- 
priate to  that  altitude.  Near  a height  h^  to  be  determined 
shortly,  the  conductivity  may  therefore  be  written 


(z-h^  /c1 

o(z,p)  =a1(h1,p)e 


(3-5) 


where  both  h^  and  are  functions  of  p. 

In  the  vicinity  of  h^,  a change  of  variable  to 

iff/ 4 (z-h1)/2^1 

y = e e 


(3-6) 
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allows  Equation  (3-4)  to  be  transformed  to 


y2  Mr  + Y + 4yoaiwci  y2^  = 0 (3_7) 

dy  y 

If  one  chooses  to  be  that  altitude  at  which 

2C2  (Poa1u))  1/2  = 1,  (3-8) 

Equation  (3-7)  takes  the  canonical  form  of  Bessel's  equation. 
The  solution  corresponding  to  outgoing  waves  at  large  alti- 
tude is 


f(p)  H 


(1) 

o 


iff/ 4 (z-h^ /2t, 

e e 


(3-9) 


Above  h-^ , which  is  the  altitude  at  which  the  skin  depth  equals 
twice  the  conductivity  scale  height,  the  potentials  fall  to 
zero  extremely  rapidly.  Thus,  h^  represents  the  upper 
boundary  of  the  waveguide,  where  ultimate  reflection  takes 
place.  It  will  be  seen  later  that  significant  reflection  also 
takes  place  at  the  altitude  hQ. 

For  altitudes  z <<  h^,  the  Hankel  function  in  Equation 
(3-9)  may  be  approximated  by  its  small  argument  expansion. 

The  scalar  potential  in  the  altitude  range  hQ  <<  z <<  h^  thus 
takes  the  form 

il»  * F(p)  [z-h1  - y-  C1]  (hQ  <<  z <<  hx)  (3-10) 

The  vector  potential  in  this  region  can  be  determined  from 
Equations  (2-17)  and  (3-1).  The  resulting  relation  is 
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iu>Az  = 


94- 

9z 


1 - 


s2(p) 


-1 


(3-11) 


Since  the  second  term  in  the  square  brackets  is  negligible  in 
this  region,  the  vector  potential  becomes 


iu)A„ 


-F(p) 


(hQ  <<  z <<  h^) 


(3-12) 


and  therefore  independent  of  altitude.  Although  it  has  not 
been  demonstrated  here,  it  can  be  shown  that  the  vector 
potential  remains  essentially  constant  over  the  entire  altitude 
range  0 <_  z <<  h^. 

It  now  remains  to  construct  a solution  which  obeys  the 
proper  boundary  condition  at  the  ground  and  has  a region  of 
validity  in  common  with  the  solution  just  derived.  As  noted 
above,  the  boundary  condition  at  the  ground  is  that  the 
horizontal  electric  field  vanishes.  The  horizontal  electric 
field  is  proportional  to  the  scalar  potential  ip , and  there- 
fore, by  virtue  of  Equation  (3-1),  to  9Az/9z.  Thus,  the 
boundary  condition  at  the  ground  is  that  9Az/3z  vanishes 
there.  With  this  boundary  condition,  a first  integral  of 
Equation  (2-17)  is 


2.  2 
n k 


A dz 
z 


(3-13) 


If  use  is  now  made  of  the  fact  that  A is  essentially  constant 

z 

in  the  range  0 < z <<  h^,  the  quantity  A may  be  taken  outside 
the  integral  for  values  of  the  upper  limit  in  this  range. 
Combined  with  Equation  (3-1),  this  leads  to 
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1 


ill 


iwAz (0 , p) 


z 


S2(p) 


dz 


(z  <<  h,  ) 
*■ 


(3-14) 


For  z >>  hQ,  the  integrand  in  Equation  (3-14)  becomes  very 
small,  and  little  error  is  made  by  extending  the  range  of 
integration  to  infinity.  Thus,  in  the  range  hQ  <<  z <<  h-^ 
the  potentials  are  given  by 


Az  = Az  (°  ' P) 


<p  = - iwAz  ( 0 , p) 


z-S  ( p) 


/ 


dz 

2 

n 


(ho  <<:  2 << 


(3-15) 


(3-16) 


The  horizontal  electric  and  magnetic  fields  derived  from 
Equation  (3-15)  and  (3-16)  must  be  made  to  agree  with  those 
derived  from  Equation  (3-12)  and  (3-13)  by  an  adjustment  of 
the  parameter  S(p).  This  is  most  easily  accomplished  by 
making  the  ratio  of  the  horizontal  electric  and  magnetic 
fields,  or  the  wave  impedance,  identical  in  the  two  cases.  It 
can  easily  be  shown  that  the  horizontal  electric  field  is  pro- 
portional to  and  the  horizontal  magnetic  field  is  proportional 
to  Az-  Thus,  the  impedance 


Z (z,p) 


Ep(z,p) 

bJTzTpT 


_ 1^  (z,  p) 

Az  ( z , p) 


(3-17) 
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and  the  requirement  that  the  impedances  be  identical  for  the 
two  solutions  results  in  the  condition 


S2(p) 


<hl  + T~  Cl> 
°°  dz 

•C  n2(z,p) 
o 


(3-18) 


For  any  given  index  of  refraction  profile,  the  eigenvalue 
S(p)  can  now  be  calculated  from  Equation  (3-18)  in  a straight- 
forward manner.  However,  it  is  possible  to  carry  the  approx- 
imate method  one  step  further.  The  integrand  in  Equation  (3-18) 
becomes  very  small  at  altitudes  very  much  above  the  height  hQ 
where 

a (h  ) = e to  ( 3-19) 

o o 


Furthermore,  the  integrand  is  essentially  unity  up  to  a 
conductivity  scale  height  or  so  below  this  altitude.  There- 
fore, an  adequate  approximation  to  the  integral  can  be 
obtained  by  writing  the  index  of  refraction  as 


= 1 + ie 


{z-ho)/Co 


(3-20) 


where  C0(p)  is  the  conductivity  scale  height  at  the  altitude 
h (p).  The  integral  in  Equation  (3-18)  can  now  be  evaluated 
analytically,  with  the  result 


(h  + i!  ; ) 

s2(p)  = — — 

<ho  - T~  to> 


(3-21) 
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The  above  is  an  approximate  analytic  expression  for  the 
eigenvalue  with  which  it  is  now  possible  to  integrate 
Equation  (2-20)  in  a rapid  and  straightforward  manner.  As 
mentioned  earlier,  the  exact  eigenvalue  S(p)  is  that  which 
would  be  obtained  for  a laterally  homogeneous  medium  with  the 
local  vertical  conductivity  profile.  A number  of  full-wave 
solutions  have  been  obtained  numerically  by  E . C.  Field  [6] 
for  a variety  of  conductivity  profiles  corresponding  to  normal 
and  disturbed  ionospheres.  In  Table  1,  the  eigenvalues 
obtained  by  Field  are  compared  with  those  calculated  from 
Equation  (3-21)  based  on  values  of  h-^,  hQ , and  £ deter- 

mined from  Field's  conductivity  profiles.  It  can  be  seen 
that  the  agreement  is  excellent  for  the  real  part  of  S(p) 
and  quite  good  for  the  imaginary  part.  In  view  of  the  basic 
uncertainties  in  conductivity  profiles,  the  use  of  the  approx- 
imation derived  above  seems  well  justified. 
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Table  1.  Calculated  Phase  Velocities  and  Attenuation 
Rates  for  a Laterally  Homogeneous  Ionosphere 


r 

v/c 

[EQ. (3-21)] 

v/c 

(EXACT) 

dB/1000  km 
[EQ.  (3-21 )] 

dB/1000  km 
(EXACT) 

10'7 

.58 

.59 

1.7 

1.8 

10'8 

.61 

.61 

1.7 

1.8 

10'9 

.65 

.65 

1.9 

1.9 

10"10 

.69 

.69 

2.0 

2.0 

10'11 

.73 

.74 

2.0 

1.95 

10"12 

.79 

.79 

1.8 

1.6 

10'13 

.81 

.805 

1.1 

1 .2 

10‘14 

.81 

.805 

0.81 

0.90 

Compari son 

of  approximate 

phase  velocities  and  attenuation  rates,  cal- 

culated  using  approximate 

Equation  (3-21),  with  Field's 

exact 

sol utions 

[6].  In  both  cases  the  conductivity  profiles 

generated  by 

Field  are 

used.  The  parameter  r is 

an  ionization  intensity  parameter 

defined  as 

r = (FY)/R2t  • 

, where  FY 

is  the  total  fission  yield  in 

megatons , 

t is  the  time  after  burst 

in  seconds,  and  R is 

the  radius  in 

kilometers 

over  which  the 

debris  is 

assumed  to  be  spread 

uni formly . 
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SECTION  IV 

SCATTERING  BY  A CYLINDRICALLY  SYMMETRIC  DISTURBANCE 


The  formalism  developed  in  the  preceding  sections  will 
now  be  applied  to  the  calculation  of  the  effect  of  a local- 
ized cylindrically-symmetric  ionospheric  disturbance  on  the 
radiated  field  of  an  ELF  transmitter.  The  applicability  of 
the  theory  requires  that  lateral  gradients  be  small;  in 
particular,  the  disturbed  region  must  join  continuously,  at 
its  outer  boundary,  with  the  undisturbed  region.  The  undis- 
turbed region  is  assumed  to  be  horizontally  stratified. 

The  geometry  of  transmitter,  disturbance,  and  receiver 
is  shown  in  plan  view  in  Figure  1.  For  simplicity,  the 
transmitter  and  receiver  have  been  assumed  to  be  located  in 
the  undisturbed  portion  of  the  waveguide.  The  coordinates 
(r,0)  are  plane  polar  coordinates  referred  to  the  transmitter 
while  (p,0)  are  plane  polar  coordinates  referred  to  the  axis 
of  the  disturbance.  The  quantities  to  be  calculated  are  the 
field  amplitudes  at  the  receiver  with  and  without  the  disturb 
ance . 

Near  the  ground,  the  vector  potential  is  a function 

g ( p , <p ) of  the  lateral  coordinates  satisfying  Equation  (2-20). 

The  function  g(p,4>)  in  the  undisturbed  region  must  satisfy 

the  same  equation  with  S = SQ.  Thus,  the  problem  may  be 

regarded  as  the  two-dimensional  scattering  by  a cylindrically 

2 2 

symmetric  complex  "potential"  V ( p ) = S (p)  - S . 

In  the  undisturbed  region,  g(p,<J>)  can  be  decomposed  into 
two  parts,  one  of  which  represents  the  incident  wave  and  the 
other  the  scattered  wave.  Thus,  g(p,4>)  can  be  written 

g(p,4>)  = g1  + gs(p,<t>)  (4-1) 
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The  function  g1  is  the  transverse  factor  of  Equation  (1-1) , with 
a replaced  by  (0-6),  i.e.. 


g1 ( r , 0 ) = cos (kScr) 


(4-2) 


By  use  of  Graf's  addition  theorem  for  Bessel  functions  [7], 

g1  can  be  expanded,  in  the  region  < p < r , in  functions 

o — — o 

of  the  ( p , 0 ) coordinates.  The  result  is 


1 ( P , <p)  = 


emJm^Sop^  t°tmcos  (m4>)  + 8msin(m<|>)] 


where 


[ p < p < r ) 

o — — o 


eo  = V2 


(4-3) 


Err,  = 1 

m 


a . . . m 
m = (-1)  cos 


(m  J 0) 


6 (kS  r ) - H (1).  (kS  r )1 

Lm+1  oo  m- 1 oo'J 


8 = (-l)m+1  sin 6/^ — )H(1)(kS  r ) 

m \ kS  r / m o o 

' o o 


(4-4) 


The  radial  function  g ( p , ) must  be  a solution  of  Bessel's 
equation  corresponding  to  outgoing  waves  at  large  radial 
distance.  The  general  outgoing  wave  solution  can  be  written 
in  the  form 


gS(p,4>)  =£3  eIT1HTi1^  (kS  P)  amcos  (m<^)  + 8 sln(m4))l 

m_  q m m o l m m mm  J 


(p  i p0> 


(4-5) 


u (1) 

where  y and 
m 

Equation  (4-3) 
region  PQ  £ P 


< 


(2) 

Y are  constants  yet  to  be  determined, 
m 

and  (4-5),  the  total  radial  function  in 
ro  is 


From 

the 


g ( p » 4>) 


X)  E (p)  cos  (m0) 
m=0  m 


+ 


( p)  sin (m0) 


] 


(4-6) 


where 


«p» 


“m[VkV>  + (kSoO'] 

Bm[Jm(kSoc)  + ’'m2>Hl1>(kSo‘»] 


(4-7) 


(pQ  £ P £ ro) 


In  order  to  complete  the  solution,  it  is  necessary  to 
solve  the  wave  equation  in  the  disturbed  region,  subject  to 
appropriate  boundary  conditions  at  p = 0 and  p = pQ . The 
function  g(p,0)  must  be  a solution  of  Equation  (2-20),  the 
general  form  of  which  is  given  by  Equation  (4-6).  The 
functions  g^^  ' ^^  ( p)  satisfy  the  equation 


1 d_ 
P dp 


(P 


+ 


k2S2(p) 


(4-8) 


(P  i PQ) 


(It  will  be  understood  that  the  omission  of  a superscript 
implies  both  superscripts.) 

The  boundary  conditions  at  p = pq  are  that  the  electric 

and  magnetic  fields  must  be  continuous  there.  It  can  easily 

2 

be  verified  that,  with  S (p)  continuous  at  p = p , all  the 
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continuity  requirements  can  be  satisfied  by  making  gm  and 
dgm/dp  continuous  at  p = pQ . 

The  boundary  condition  at  p = 0 is  that  all  the  fields 
remain  finite,  which  in  turn  requires  9m(0)  and  g^(0)  to  be 
finite.  This  boundary  condition  can  be  conveniently  formal- 
ized by  transforming  Equation  (4-8)  into  a Ricatti-type 
equation.  First  one  lets 

gm(p)  = pm  Um(p)  (4-9) 


where 


P = P/P0 

The  function  U(p)  then  satisfies 


u"  + ( 2m+l)  Ug,  + £2  2 = 

m p m 


(4-10) 


where 


k = kp. 


If  one  then  introduces  as  the  variable 


Y — _ s(2m+l)  jti 
m M U 

m 


Equation  (4-10)  becomes 


' r2„2 - (2m+l)  , m 

Ym  = k S (p)p  + ZTJ^T) 


(4-11) 


(4-12) 
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This  equation  can  be  integrated  outward  from  the  origin  in 

a straightforward  manner  once  the  behavior  of  Y at  the 

m 

origin  is  known.  This  can  be  deduced  from  the  requirement 
that  gm(0)  and  9m(0)  are  finite.  In  the  neighborhood  of 
P = 0,  S2(p)  is  approximately  constant,  and  Equation  (4-8) 
reduces  to  Bessel's  equation.  The  solution  which  is  finite 
at  the  origin  varies  with  p as 

gm  s:  pm[l  + 0 ( p2 ) ] p ■+  0 (4-13) 

Thus,  from  Equation  (4-9)  and  (4-11) 

Ym(P)  = P2  (m+1)  P ->  0 (4-14) 


so  that 


Ym(0)  = Y>)  = 0 


(4-15) 


The  last  condition  is  sufficient  to  start  an  outward  numerical 
integration  of  Equation  (4-12)  , from  which  one  obtains  a value 

Y ( i_ ) at  the  boundary  p = p . 

m ° ( 1)  12) 

Finally,  the  constants  y^  and  y^  in  Equation  (4-7) 

can  now  be  determined  from  the  boundary  conditions  at  p = pQ. 

Since  9m(p)  and  g^(p)  must  be  continuous  at  p = pQ,  their 

ratio  must  also  be  continuous.  In  the  undisturbed  region, 

the  ratio  is  calculated  from  Equation  (4-7),  and  in  the 

disturbed  region  from  Equation  (4-9)  and  (4-11).  The  resulting 

relation,  for  either  superscript,  is 


kSo 


J ' (kS  p ) + Y H ( 1)  (kS  p ) 
L m o o mm  o o . 

J (kS  p ) + Y Ili,1*  (kS  o )1 
m oo  mm  ool 


m - Y (1) 
m 


(4-16) 
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r ^ 


The  constants  v ^ and  y ^ are  clearly  the  same,  and  with  a 
mm 

little  further  manipulation,  can  be  shown  to  be  given  by 

Ym  " [kSo(,oHl"(kSof'o>-V1»H'!'1,,kSopo)] 

Combined  with  Equation  (4-5) , this  gives 

oo 

gS  ( p , (p)  = ^ emYmH^1^  (kSQp)  [amcos  (m<})) +3msin  (m<{))  ] 
m=  0 

which  now  enables  one  to  calculate  the  scattered  amplitude 
everywhere . 


(4-17) 


(4-18) 


SECTION  V 


CALCULATION  OF  SCATTERED  MAGNETIC  FIELD 


In  ELF  reception  by  receivers  submerged  in  the  ocean,  the 
receiving  antenna  measures  the  horizontal  electric  field.  The 
horizontal  electric  field  at  any  depth  in  the  ocean  can  be 
determined  from  the  horizontal  magnetic  field  at  the  surface. 

It  is  therefore  the  latter  quantity  which  will  be  calculated 
explicitly  here. 

It  can  be  expected  that  an  ionospheric  disturbance  will 
have  its  greatest  effect  when  its  axis  passes  through  the 
transmitter-receiver  line.  Consequently,  it  will  be  assumed 
that  0 = 0 and  <J>  = 0 , tt  (see  Figure  1)  . Furthermore,  it  will  be 
assumed  that  the  transmitter  is  oriented  for  maximum  signal 
at  the  receiver,  so  that  6=0.  These  assumptions,  in  addi- 
tion to  applying  to  the  case  of  greatest  practical  interest, 
also  simplify  the  calculation  somewhat.  However,  the  general- 
ization to  arbitrary  configuration  of  transmitter,  receiver, 
and  disturbance  is  straightforward. 

For  the  configuration  described  above,  the  horizontal 
magnetic  field  can  be  shown  to  have  only  an  azimuthal  component. 
The  azimuthal  magnetic  field  is  related  to  the  vector  potential 
by 

SA 

B,n  (5-1) 

Aside  from  an  excitation 
field  is  therefore 


factor,  the  total  azimuthal  magnetic 


E (-i 


m=0 


3g 

m ^m 

e a t; 

m m 9p 


(5-2) 


(r  < r0> 
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m=0 


9g 


1 mam  9p 


m 


(5-3) 


(r  1 ro} 


where  a is  given  by  Equation  (4-4)  with  6=0.  (The  factor 
m ^ 

(-l)n  which  appears  in  Equation  (5-2)  arises  from  the  factor 


cos(ni<(>)  in  Equation  (4-6)  , in  which  <f>  = ti  for  r < rQ.) 
azimuthal  field  of  the  incident  wave  is  given  by 


The 


~ H.(1)  (kS  r) 
9p  1 o 


(5-4) 


r = ro  T p 


with  the  upper  sign  corresponding  to  Equation  (5-2)  and  the 


lower  sign  to  Equation  (5-3). 


The  evaluation  of 


m 


9 P 


in  Equation  (5-2)  and  (5-3)  depends 


on  whether  the  receiver  is  inside  or  outside  the  disturbed 
region.  If  the  receiver  is  outside  the  disturbed  region 
(P  > PQ)  , 


gm(P)  = (kS  P)  + J (kSp) 


(5-5) 


where  y is  calculated  from  Equation  (4-17)  after  an  integration 
m 9g 

of  Equation  (4-12)  from  p = 0 to  p = 1 . The  quantity  — can 

then  be  calculated  by  a straightforward  differentiation  of 

Equation  (5-5) . 

If  the  receiver  is  inside  the  disturbance  (p  _<  p ) , one 

3gm 

additional  integration  is  required  to  determine  -yj-.  It  is 
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easy  to  show  from  Equation  (4-9)  and  (4-11)  that  gm  satisfies 
the  differential  equation 


g,>>  - gm(5) 


m + V5’ 


. 2m+l 


(5-6) 


with  the  boundary  condition  (Equation  (5-5) 


g (1)  = YmH (1)  (kS  ) + J (kS  ) 
m m o m o 


(5-7) 


The  function  Ym(P)  can  be  tabulated  during  the  same  integration 
of  Equation  (4-12)  from  which  y is  calculated.  With  y and 
Ym(p)  known,  Equation  (5-6)  can  be  integrated  inward  from 


p = 1 to  determine 


’m 


9p 


anywhere  inside  the  disturbance. 


For  a configuration  in  which  the  transmitter  is  under  the 
center  of  the  disturbance  and  the  receiver  is  outside  the  dis- 
turbance, the  field  at  the  receiver  can  be  obtained  by  inter- 
changing the  positions  of  transmitter  and  receiver,  performing 
the  calculation  as  described  above,  and  invoking  reciprocity. 
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SECTION  VI 


NUMERICAL  RESULTS  AND  DISCUSSION 


Illustrative  calculations  based  on  the  preceding  theory 
have  been  carried  out  for  the  disturbances  resulting  from  a 
high-altitude,  megaton-range  nuclear  explosion  during  daytime. 
The  bomb  fission  yield  was  taken  to  be  2 MT  and  the  fission 
products  were  assumed  to  rise  to  an  altitude  of  300  km  at 
1 min  and  1000  km  at  10  min  after  the  burst.  For  both  times, 
the  ionization  rate  from  delayed  gamma  rays  was  calculated  as 
a function  of  altitude  and  distance  from  ground  zero  using 
standard  formulas  [8] . Ionization  from  delayed  beta  rays  was 
neglected,  since  it  depends  on  the  local  magnetic  dip  angle 
and  is  expected  to  have  a smaller  effect  than  the  gammas. 

The  resulting  electron  and  ion  concentrations  were 
calculated  by  R.  Turco  [9]  assuming  steady-state  chemistry 
and  including  e , NO+.nH2<0,  H^O+.nH20,  02  and  NO^.nH-O.  Using 
mass-dependent  ion  mobilities  [10]  similar  to  those  of  Carroll 
and  Mason  [11]  , ionospheric  conductivity  profiles  were 
computed  [9].  From  these  results,  the  two  altitudes  and  the 
two  scale  heights  needed  for  the  calculation  of  approximate 
eigenvalues  (Equation  (3-21))  at  45  Hz  were  determined.  The 
results  appear  in  Tables  2 and  3,  the  last  entry  in  each  table 
representing  the  beginning  of  ambient  conditions.  These  tables 
were  the  input  data  for  electromagnetic  field  calculations 
using  the  methods  described  in  the  preceding  section. 

In  each  case,  the  transmitter  was  located  outside  the 
disturbance  at  1000  km  from  the  nearest  edge,  and  the  electro- 
magnetic field  was  calculated  at  points  along  the  line  between 
the  transmitter  and  the  center  of  the  disturbance,  ranging 
from  200  km  from  the  transmitter  to  1000  km  beyond  the  far  edge 
of  the  disturbance.  (Larger  transmitter-receiver  distances 
would  require  earth-curvature  corrections.)  The  fields 
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Table  2.  Ionospheric  Parameters  and  Approximate  Eigenvalues  for  2 MT 
of  Fission  Products  at  300  km  Altitude  1 Min  After  Burst  Time 


MiMlIJIi! 

10021 

m 

Re(S) 

Im(S) 

0 

18.66 

3.19 

59.32 

2.97 

1 .73 

.30 

200 

19.60 

3.18 

60.22 

2.92 

1.71 

.28 

400 

22.60 

3.54 

62.47 

2.73 

1.62 

.25 

600 

25.19 

3.44 

63.88 

2.78 

1.56 

.22 

800 

27.76 

3.71 

65.03 

2.79 

1.50 

.21 

1000 

29.88 

3.47 

65.95 

2.77 

1.46 

.13 

1200 

31.87 

3.42 

66.71 

2.74 

1.43 

.17 

1400 

34.30 

3.23 

67.34 

2.71 

1.39 

.15 

1600 

36.70 

3.01 

67.62 

2.69 

1.35 

.13 

1800 

40.10 

2.64 

68.04 

2.76 

1.30 

.11 

2000 

44.58 

2.19 

68.50 

2.77 

1.24 

.087 

2200 

50.86 

1.73 

69.85 

2.72 

1.17 

.067 

2300 

54.63 

1.65 

71.17 

2.78 

1.14 

.062 

2400 

58.13 

2.11 

73.68 

2.92 

1.12 

.067 

2450 

58.28 

2.50 

75.53 

2.93 

1.14 

.073 

2500 

58.20 

2.57 

77.58 

3.29 

1.15 

.078 

2550 

58.20 

2.57 

78.24 

3.42 

1.16 

.080 

2600 

58.20 

2.57 

78.20 

3.52 

1.16 

.081 

Table  3.  Ionospheric  Parameters  and  Approximate  Eigenvalues  for  2 MT 
of  Fission  Products  at  1000  km  Altitude  10  Min  After  Burst  Time 


GROUND 

RANGE 

(km) 

ho 

(km) 

cO 

(km) 

hl 

(km) 

C1 

(km) 

Re(S) 

Im(S) 

0 

33.70 

8.52 

70.05 

3.15 

1.36 

.31 

200 

33.81 

8.49 

70.08 

3.16 

1.36 

.31 

400 

34.45 

8.31 

70.27 

3.21 

1.35 

.30 

700 

35.80 

8.14 

70.67 

3.33 

1.34 

.28 

1000 

37.29 

8.39 

71.11 

3.47 

1.31 

.28 

1500 

40.00 

7.15 

71.77 

3.71 

1.30 

.23 

2000 

41.88 

6.66 

72.38 

4.01 

1.28 

.21 

2600 

44.13 

5.34 

73.10 

4.12 

1.27 

.18 

3200 

47.13 

3.90 

73.72 

4.20 

1.24 

.14 

3400 

49.12 

3.12 

73.90 

4.21 

1.22 

.12 

3500 

50.44 

2.68 

73.99 

4.22 

1.21 

.10 

3600 

52.23 

2.23 

74.12 

4.20 

1.19 

.094 

3800 

57.21 

1 .87 

75.69 

3.86 

1.15 

.076 

4400 

58.20 

2.57 

78.19 

3.52 

1.16 

.081 
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corresponding  to  the  first  ten  partial  waves  were  calculated, 
and  superimposed  to  give  the  scattered  field  at  any  point. 

In  all  cases,  only  the  lowest  three  or  four  partial  waves 
were  found  to  make  an  appreciable  contribution  at  any  field 
point . 

In  Figures  2 and  3 the  attenuation  of  the  45-Hz  elec- 
tromagnetic field  is  plotted  as  a function  of  distance  from 
the  transmitter.  The  attenuation  in  the  absence  of  the 
disturbance  is  also  plotted  for  purposes  of  comparison  at 
1 min  after  the  explosion  (Figure  2) ; the  radius  of  the  dis- 
turbance is  about  2600  km,  and  the  ionosphere  at  the  center 
has  been  lowered  by  40  km.  The  additional  attenuation  across 
the  entire  disturbance  is  about  0.5  dB.  At  10  min  after  the 
explosion  (Figure  3)  the  radius  of  the  disturbance  is  about 
4400  km,  the  ionosphere  at  the  center  has  been  lowered  by 
25  km,  and  the  additional  attenuation  is  about  3 dB.  In  both 
cases,  the  field  between  the  transmitter  and  the  leading  edge 
of  the  disturbance  is  essentially  unaltered,  indicating  that 
there  is  little  back-scattering. 

It  is  interesting  to  compare  these  results  of  the  full- 
wave  solution  with  the  "great-circle  WKB"  approximation  which 
has  been  used  in  ELF  calculations  by  Pappert  and  Moler  [12] 
and  others.  This  is  essentially  a two-dimensional  approxima- 
tion in  which  the  properties  of  the  waveguide  are  allowed  to 
vary  along  the  direction  of  propagation,  but  are  assumed  uniform 
in  the  perpendicular  direction.  The  properties  throughout  the 
waveguide  are  assumed  to  be  those  of  the  great-circle  path 
between  transmitter  and  receiver.  In  this  approximation,  the 
magnetic  field  of  the  wave  is  given  by 


B , ( r ) 


.1/2 


[hoVr)] 


r)l1/2 


T 

S1  S (r) 


-1/4 


exp 


[-/ 


Im(s)dr 


(6-1) 
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where  the  superscript  T denotes  ti.e  transmitter  location. 

-1/2 

The  factor  [h  (r) ] represents  the  increase  or  decrease 

in  energy  density  as  the  height  of  the  waveguide  changes. 

The  attenuation  calculated  from  this  approximation  is 
also  shown  in  Figures  2 and  3.  The  approximation  has  been 
applied  only  to  that  portion  of  the  path  inside  the  disturbance, 
so  that  Equation  (6-1)  has  been  normalized  to  the  undisturbed 
value  at  the  leading  edge.  The  results  show  that  this  approx- 
imation underestimates  the  attenuation  in  most  of  the  disturb- 
ance region,  but  overestimates  the  attenuation  at  the  far  edge 
of  the  disturbance  by  about  a factor  of  two.  This  behavior 
is  not  surprising,  since  the  two-dimensional  approximation 
forces  the  wave  to  go  under  the  ionospheric  depression,  thus 
both  compressing  and  absorbing  it,  while  the  realistic  full- 
wave  solution  permits  the  radio  wave  to  go  around  the 
depression . 
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